Effective-Hamiltonian modeling of external pressures in ferroelectric perovskites 
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The phase-transition sequence of a ferroelectric perovskite such as BaTiOa can be simulated 
by computing the statistical mechanics of a first-principles derived effective Hamiltonian [Zhong, 
Vanderbilt and Rabe, Phys. Rev. Lett. 73, 1861 (1994)]. Within this method, the effect of an 
external pressure (in general, of any external field) can be studied by considering the appropriate 
"enthalpy" instead of the effective Hamiltonian itself. The legitimacy of this approach relies on two 
critical assumptions that, to the best of our knowledge, have not been adequately discussed in the 
literature to date: (i) that the zero-pressure relevant degrees of freedom are still the only relevant 
degrees of freedom at finite pressures, and (ii) that the truncation of the Taylor expansion of the 
energy considered in the effective Hamiltonian remains a good approximation at finite pressures. 
Here we address these issues in detail and present illustrative first-principles results for BaTiC>3. 
We also discuss how to construct effective Hamiltonians in cases in which these assumptions do not 
hold. 
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I. INTRODUCTION 

A few years ago it was shown that it is possible to re- 
produce the phase transition sequence of a ferroelectric 
such as BaTi03 by simulating the statistical mechanical 
properties of an effective Hamiltonian built on the ba- 
sis of first-principles calculations p]. Since then, this 
approach has been used to study the phase diagrams 
and electromechanical responses of other perovskite crys- 
tals (SrTi0 3 @, PbTi0 3 [§, and KNb0 3 {§), and even 
of disordered perovskite solid solutions (Pb(Zri_ x Ti a; )03 
(PZT) g and Pb(Sc .5Nbo.5)0 3 @). These calculations 
have typically resulted in a very good qualitative agree- 
ment with experiment and clear physical pictures of the 
phenomena studied. Presently, even though the method 
usually gives poor results for the actual values of the tran- 
sition temperatures and other quantities, it is regarded as 
having the predictive power necessary to tackle materials- 
design problems 0, || . 

An effective Hamiltonian can be viewed as a Taylor 
expansion of the energy of the system, in terms of a set 
of relevant degrees of freedom, around a high-symmetry 
reference structure. For ferroelectric perovskites, the rel- 
evant variables are typically the strains and, most im- 
portantly, the local polar distortions that sum up to pro- 
duce the spontaneous polarization. The high-symmetry 
reference structure is chosen to be the non-polar, cubic 
perovskite phase, whose equilibrium volume is computed 
ab initio. 

Very conveniently, within this approach it is straight- 
forward to consider the effect of external fields on the sys- 
tem. One just has to add a term to the Hamiltonian that 
transforms the energy into the appropriate enthalpy. For 
instance, an external hydrostatic pressure p would enter 
the model in a term pV, where V is the volume of the 
simulated system. Examples of such applications are the 
calculation of the T-p phase diagram of BaTiC>3 || and, 
very recently, the study of electric-field driven transition 
paths in PZT @. 



This approach to modeling external fields stems di- 
rectly from statistical mechanics, and it is exact provided 
the considered Hamiltonian is the complete Hamiltonian 
of the system. Using effective incomplete Hamiltonians 
actually involves a number of implicit approximations 
in the description of the effect of an external field. To 
the best of our knowledge, these implicit approximations 
have neither been adequately discussed, nor has their im- 
portance been carefully tested, in the literature to date. 

In this paper we address these questions in detail. For 
concreteness, we consider the case of an external pressure, 
but our formal argument applies to any external field. We 
illustrate the discussion with calculations of the effect of 
pressure on BaTiC>3, which should be representative of 
ferroelectric perovskites and is specially interesting be- 
cause of the discrepancy between the theoretical || and 
experimental jllj] results for its T-p phase diagram. 

The paper is organized as follows. In Section II we dis- 
cuss the pressure dependence of the relevant degrees of 
freedom that define the effective Hamiltonian. In Sec- 
tion III we show how the truncation in the effective- 
Hamiltonian energy expansion actually implies an ap- 
proximation in what we call the "pressure dependence" 
of the parameters of the model. In Section IV we tackle 
the issue of the anharmonic couplings between the cho- 
sen relevant degrees of freedom and the rest of variables 
in the system. Finally, we summarize and present our 
conclusions in Section V. 



II. RELEVANT DEGREES OF FREEDOM AS A 
FUNCTION OF PRESSURE 

The fundamental variables entering effective Hamilto- 
nians for ferroelectric perovskites are the localized, polar 
displacement patterns associated with the spontaneous 
polarization. In order to define these local polar modes 
ab initio, one examines the calculated phonon [ jlll dis- 
persion curves of the high-symmetry cubic phase; the soft 
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phonons (i.e., those with very small or negative force con- 
stant) are obviously those to be included in the model. 
(Note that the spontaneous polarization will correspond 
to an unstable zone-center phonon.) The local modes, 
or lattice Wannier functions, are then obtained from the 
corresponding eigenvectors Jl3| |. 

Taking the calculated equilibrium, high-symmetry 
phase as a reference implies that the relevant variables we 
determine are actually the relevant variables at zero pres- 
sure. Including the effect of pressure by simply adding 
an extra pV term to the Hamiltonian thus relies on the 
fundamental assumption that the relevant variables at 
zero pressure will remain the relevant variables at finite 
pressure. There are two situations in which this assump- 
tion could fail, namely the applied pressure could give 
raise to additional soft modes and/or the zero-pressure 
relevant configuration space could change significantly. 
In the former case, the obvious solution is to include the 
new soft modes in the model. In the latter, one would 
have to redetermine the relevant polar degrees of free- 
dom, now under applied pressure, and recalculate the 
effective Hamiltonian accordingly. 

We have performed first-principles calculations to 
study these issues in BaTiC>3 Q. In this material, the 
modes that most likely may become soft under compres- 
sion are those at the M and R zone-boundary points that 
involve rotations of the oxygen octahedra. We have cal- 
culated the evolution of the eigenvalues of such modes 
as a function of the lattice parameter, ranging from a = 
7.46 a. u. (the calculated equilibrium lattice parameter of 
the cubic phase) to a = 7.36 a. u. (which is well beyond 
the value at which the ferroelectric instability disappears; 
see next section). We find no softening with increas- 
ing pressure; on the contrary, the modes become slightly 
harder. (For instance, the eigenvalues of the modes at 
R, which are softer than those at M, are 0.0200 a. u. and 
0.0220 a. u. for a = 7.46 a. u. and a = 7.36 a. u. respec- 
tively.) On the other hand, the evolution with pressure 
of the relevant variables related to the ferroelectric insta- 
bility can be monitored by observing the change of the 
zone-center soft mode. We find that the overlap between 
the unstable eigenmode at a = 7.46 a.u. and the corre- 
sponding, nearly-unstable eigenmode at a = 7.36 a.u. is 
above 98%. In summary, for BaTiC>3 the relevant vari- 
ables defined at the zero-pressure equilibrium lattice con- 
stant remain a very good approximation all through the 
interesting pressure range. 



III. "PRESSURE DEPENDENCE" OF THE 
HAMILTONIAN PARAMETERS 

Let us assume that the zero-pressure relevant degrees 
of freedom continue to be the relevant ones in the pres- 
sure range of interest. In the following we argue that, 
even in this case, the effect of pressure is not fully cap- 



tured by the addition of the pV term. 

For simplicity, we consider a one dimensional ferroelec- 
tric, and describe it in terms of the macroscopic variables 
P (polarization) and 77 (strain). For small P and 77, the 
energy of the system can be written as a low-order Taylor 
expansion around the high-symmetry paraelectric phase 
(P = r/ = 0), 

E{P, T])^Eo- C 2 . P 2 + Ci,oP 4 + C ,2V 2 ~ C 2 .iP 2 V, 

(1) 

where all the expansion coefficients (in obvious notation) 
are chosen to be positive. (The maximum orders included 
in the expansion of Eq. (|l]) are those normally adopted in 
the effective Hamiltonians proposed to date in the litera- 
ture ||, ||, |(| . ) Let us include the effect of an external 
pressure p by adding the term pVo(l + rj), with Vq the 
equilibrium volume of the paraelectric phase. It is now 
convenient to consider the change of variables 
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-p = 77 + rj p 



(2) 



where Vb(l + rj p ) is the equilibrium volume of the para- 
electric phase under applied pressure. The enthalpy then 
becomes 



E(P, mP ) = E (p) 



„ VoC2 } l 

C2.0 - — p 



P z 



-CaqP 



C ,2^ - C 2A P 2 ?j, (3) 



where Eq (p) is the pressure-dependent energy of the non- 
polar phase. If we now denote the bracketed prefactor 
of the P 2 term as 6*2,0, and also introduce Eq, we can 
rewrite the enthalpy as 



E(P, ?j)=E - C 2 .qP 2 + C ifi P 4 + Co,2Tp - C 2 . X P% : 



(4) 



which is formally identical to Eq. (|I|). Thus, we have 
proven that, within our model, the system under ap- 
plied pressure is in all respects equivalent to the origi- 
nal, pressure-free system described by Eq. (|l|), but with 
some modified parameters (indicated by tildes) that can 
be regarded as being "pressure dependent" . (In the fol- 
lowing, references to the "pressure dependence" of the 
parameters of the Hamiltonian should be understood in 
this sense.) 

Having reformulated the problem in this way, an obvi- 
ous question arises. If Eq and C 2t Q are modified by pres- 
sure, why is it that the remaining parameters in Eq. (|l|) 
are not? The answer is that pressure should indeed affect 
the rest of the parameters in Eq. (|l|) , but this effect is not 
captured by our model because we are using a truncated 
Taylor expansion of the energy. Imagine, for instance, 
that Eq. (|l|) contained the higher-order term C^.\P^r\. 
The variable change in Eq. (0) would then generate new 
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FIG. 1: Calculated eigenvalue of the ferroelectric soft mode of 
BaTiOs as a function of pressure (dots), as compared with the 
behavior predicted (see text) from the zero-pressure effective 
Hamiltonian (solid line). (See E3| for the technical details 
of the ab initio calculations.) From the explicit calculation 
we obtain a ferroelectric transition point at p c — 7.0 GPa 
(Vc/Vo — 0.966), while the effective-Hamiltonian approach 
places the transition at p c — 6.7 GPa (V c /Vb = 0.964, as given 
by the equation of state built into the effective Hamiltonian 
itself). In both cases, the calculated equilibrium volume of 
the cubic phase is Vo = 414 a.u. 

terms CuP^rj + C4 l^pP 4 , thus modifying 64,0 to be- 
come C\o = C4 j0 + C^iTJp. 

The previous argument can readily be applied to the 
microscopic effective Hamiltonians of interest ]l7[| . We 
have thus shown that a truncation of the energy expan- 
sion, which may be perfectly well justified for the cal- 
culation of zero-pressure properties, actually implies an 
approximation in the pressure dependence of the param- 
eters of the model. Obviously, we can improve on this 
approximation by including higher-order terms in the ef- 
fective Hamiltonian. However, our experience tells us 
that it is usually difficult to calculate high-order terms 
accurately, and the very issue of where to truncate the ex- 
pansion may become a problem. There is, fortunately, an 
alternative approach suggested by the derivation above. 
One can instead take the equilibrium, high-symmetry 
phase at pressure p as the reference structure from which 
one calculates the Taylor expansion of the energy. Such 
an expansion, which could likely be kept low-order, would 
be analogous to Eq. ([|) but with all the parameters con- 
taining their full pressure dependence by construction. 

In order to make these ideas explicit, we have calcu- 
lated the pressure dependence of the parameters in the 
effective Hamiltonian of BaTiOs juj. The soft-mode 
eigenvalue, which is exactly the equivalent of — 2C2,o in 
Eq. (Q), is computed as a function of the unit-cell vol- 
ume, and the corresponding pressure is determined from 
the equation of state of the cubic phase. The results are 



shown in Fig. [|. We find that the linear law resulting 
from the pV term (see the bracketed term in Eq. (0)) 
turns out to agree surprisingly well with the values ex- 
plicitly calculated at different pressures. We also com- 
puted the pressure dependence of the parameters that are 
analogous to 6*4,0, Co, 2 and 62,1 in Eq. (|l]). In all cases 
the values calculated as a function of pressure lay within 
approximately 10% of the zero-pressure value. Very im- 
portantly, this applies to the parameter that determines 
the relative stability of the different ferroelectric phases 
of BaTi0 3 (7' in Ref. @). Therefore, we find that the 
ground state of the system is rhombohedral throughout 
the pressure range in which the ferroelectric instability 
exists, in agreement with the first-principles effective- 
Hamiltonian results of Zhong et al. B . 

IV. THE COUPLING BETWEEN RELEVANT 
AND IRRELEVANT VARIABLES 

The experimental work of Ref. [jllj suggests that, along 
isotherms at very low temperature, BaTiOs undergoes 
a transition sequence with increasing pressure that pro- 
gresses from the zero-pressure ferroelectric rhombohedral 
phase to ferroelectric orthorhombic, then ferroelectric 
tetragonal, and finally paraelectric cubic phases. As we 
have seen in the previous sections, we have found no in- 
dication of such a behavior in our first-principles results. 
This disagreement is not necessarily a serious one, since 
the authors of Ref. |llj] attributed their observed transi- 
tion sequence to quantum fluctuations (zero-point atomic 
motion), an effect that is not taken into account in our 
calculations. In any case, it led us to explore yet an- 
other possible effect of pressure that could, in principle, 
account for such a discrepancy. 

The effective-Hamiltonian approach relies on the as- 
sumption that the coupling between the chosen relevant 
variables and the rest of the degrees of freedom in the 
system is not important. Since the relevant modes are de- 
fined in terms of phonons that diagonalize the harmonic 
part of the total Hamiltonian of the system, such a cou- 
pling is anharmonic, which partly justifies our neglecting 
it. However, we know that in BaTiOs, a s we approach 
the critical pressure at which the ferroelectric instabil- 
ity disappears, the energy differences among the various 
low-symmetry phases go to zero. Small anharmonic ef- 
fects could then be essential to the determination of the 
ground state of the system. Note that, if this were the 
case, the usual effective-Hamiltonian approach, in which 
the relevant variables are identified by looking at the har- 
monic phonon dispersion curves, should be modified. 

We can study this issue very easily in BaTiOs by fully 
relaxing the structure of the system in different con- 
strained symmetries, while keeping the volume fixed [ plj| . 
The symmetries we consider are of course the tetragonal, 
orthorhombic, and rhombohedral ones corresponding to 
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FIG. 2: Energies of the different phases of BaTiOs (cubic, 
tetragonal, orthorhombic, and rhombohedral) as a function 
of the volume of the unit cell. The energies are computed by 
relaxing the atomic positions and cell shape under the con- 
straints of preserving the corresponding symmetry and keep- 
ing the cell volume fixed. (See [u| for the technical details.) 
We have marked with a vertica line the ferroelectric transi- 
tion point between the non-polar cubic and polar rhombohe- 
dral phases, located at V c /Vo = 0.974 (p c = 5.6 GPa) where 
Vo = 417 a.u. is the calculated equilibrium volume of the cubic 
phase. The discrepancy between these numbers and those in 
the caption of Fig. ^ is not surprising if we take into account 
that they have been computed using two different codes (and 
two different kinds of pseudopotential schemes). 



the ferroelectric phases of the material, and we explore 
the interesting range of volumes. Figure ^ shows the en- 
ergy of the relaxed structures as a function of volume. 
For V/Vq < 0.980 the energies of the different phases dif- 
fer by less than 0.01 mHa. Such small differences are be- 
yond the accuracy of our computational technique and, 
therefore, we are not able to study the critical region 
in detail. Nevertheless, we can determine the critical 
volume quite reliably as the point at which 6*2,0 passes 
through zero; we obtain V c /V~o = 0.974, indicated with 
a vertical line in Fig. ||. On the other hand, we find 
that, down to V/Vq ~ 0.977, there is no common tan- 
gent among the curves in Fig. ^ that would correspond 
to a first-order phase transition between two ferroelectric 
phases at p < p c . It thus seems reasonable to conclude 
that the calculated ground state is rhombohedral up to 
p c , the point at which it undergoes a second-order tran- 
sition to the paraelectric cubic phase. This result is in 
complete agreement with the predictions of the effective- 
Hamiltonian approach. 

The importance of the above-mentioned anharmonic 
couplings in determining the ground state of the sys- 
tem can be quantified, for a given volume, by expand- 
ing the displacement vector corresponding to the relaxed 



atomic positions as a linear combination of the zone- 
center modes calculated at that same volume. We find 
that, throughout the interesting range, the ferroelectric 
soft mode accounts for approximately 99% of the total 
atomic relaxation. This result further justifies neglecting 
the anharmonic couplings among relevant and irrelevant 
variables in the case of BaTiC>3. 



V. SUMMARY AND CONCLUSIONS 

In the usual approach to simulating external pressures 
within the effective-Hamiltonian scheme, one starts from 
the calculated zero-pressure effective Hamiltonian and 
adds to it a pV term to obtain the appropriate enthalpy. 
Such a procedure implies a number of approximations 
that we have discussed in detail in this paper. Our dis- 
cussion is actually general, and applies to any external 
field. 

On the one hand, we have seen that there are approx- 
imations related with the determination of the relevant 
degrees of freedom to be included in the model. In princi- 
ple, one should determine the relevant variables as a func- 
tion of pressure. One should also check that the relevant 
variables do not have significant anharmonic couplings 
with the rest of the degrees of freedom of the system 
throughout the interesting pressure range. 

On the other hand, we have discussed a rather subtle 
implicit approximation that has to do with the trunca- 
tion of the energy expansion in the effective Hamiltonian. 
We have shown that such a truncation, which may be jus- 
tified for the purpose of zero-pressure calculations, actu- 
ally implies an approximation in what we have called the 
"pressure dependence" of the parameters of the model. 
As a solution to this problem we propose the calculation 
of effective Hamiltonians "under pressure," meaning that 
for a given pressure we would take as our reference struc- 
ture the calculated, equilibrium high-symmetry phase at 
that pressure. 

We have studied in detail all these issues for BaTi03. 
We have shown that all the previous approximations are 
actually highly accurate for this material, which we be- 
lieve should be representative of ferroelectric perovskites. 
The usual approach to simulating external pressures thus 
seems to be fairly reliable. 

Regarding the T-p phase diagram of BaTiOs, our re- 
sults indicate that at a classicallevel (i.e., without consid- 
ering the zero-point motion of the ions) the only possible 
polar phase at K is the rhombohedral one. On the other 
hand, as Fig. ^| shows, the differences in energy among 
the ferroelectric phases become very small as the volume 
decreases towards V c . It thus seems possible that, at low 
temperatures, the quantum fluctuations of the ions may 
be able to change the relative stability of the ferroelectric 
phases and lead to a phase diagram similar to the one re- 
ported by Ishidate et al. (see Fig. 1 in Ref. 111]]). This 
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very interesting question will be addressed elsewhere . 

We thank Alberto Garcia and Laurent Bellaiche for 
their help regarding various technical issues. This work 
was supported by the ONR through Grant N00 14-97-1- 
0048. 



13236 (1995); and in J. 
Perez-Mato, Phys. Rev. 



Iniguez, A. 
B 61, 3127 



1] W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. 

Lett. 73, 1861 (1994). 
2] W. Zhong and D. Vanderbilt, Phys. Rev. Lett. 74, 2587 
(1995). 

3] U. V. Waghmare and K. M. Rabe, Phys. Rev. B 55, 6161 
(1997). 

4] H. Krakauer et al, J. Phys. Cond. Matter 11, 3779 
(1999). 

5] L. Bellaiche, A. Garcia, and D. Vanderbilt, Phys. Rev. 

Lett. 84, 5427 (2000). 
5] R. Hemphill, L. Bellaiche, A. Garcia, and D. Vanderbilt, 

Appl. Phys. Lett. 77, 3642 (2000). 
7] J. Iniguez and L. Bellaiche, Phys. Rev. Lett. 87, 095503 
(2001). 

8] A. M. George, J. Iniguez, and L. Bellaiche, Nature 413, 
54 (2001). 

9] W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. 

B 52, 6301 (1995). 
0] L. Bellaiche, A. Garcia, and D. Vanderbilt, Phys. Rev. B 

64, 060103 (2001). 
1] T. Ishidate et al, Phys. Rev. Lett. 78, 2397 (1997). 
2] All through the paper we use the terms phonon or normal 
mode to refer to the eigenvectors of the force- const ant 
matrix. 

[13] Different ways of calculating the local modes can be 
found in Ref. M; in K. M. Rabe and U. V. Waghmare, 



[1] 
[2] 
[3] 
[4] 
[5] 
[6] 
[7] 
[8] 
[9] 
[10] 

[11] 
[12] 



[14] 



[15] 
[16] 

[17] 



[18] 



Phys. Rev. B 52, 
Garcia, and J. M. 
(2000). 

For these calculations we used density-functional the- 
ory within the pseudopotential approach (using Vander- 
bilt 's ultrasoft pseudopotentials |15j|) and the local den- 
sity approximation. We used a 25 Rydberg cutoff for 
the plane-wave basis and a 6 x 6 x 6 Monkhorst-Pack 
special k-point mesh. Importantly, in our calculations 
we considered 10 valence electrons for Ba (5s 2 5p 6 6s 2 ), 
12 for Ti (3s 2 3p 6 3d 2 4s 2 ), and 6 for O (2s 2 2p 4 ). These 
are well- converged calculation conditions, as discussed in 
Ref. Jl|]. 

D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 

R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 49, 

5828 (1994). 

In the typical effective Hamiltonian, where the strain 
is three-dimensional, the pV factor actually produces a 
term analogous to Co, 2 as well as other strain-related 
changes. However, the high-order terms associated with 
the polar modes, and the terms that couple them with 
the strains, are not changed. 

We performed these calculations using the Vienna Ab- 



[19] 

[20] 
[21] 



19| 



We used a 6 x 
aerg cutoff for the 



initio Simulation Package (VASP) 
6x6 special k-point grid and a 36 Ryd 
plane-wave basis, and we worked within the local density 
approximation. We used the Projector Augmented Wave 
(PAW) method @, and valence electrons as indicated 
in @. 

G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993); 
G. Kresse and J. Furthmiiller, Comput. Mat. Sci. 6, 15 
(1996). 

P. Blochl, Phys. Rev. B 50, 17953 (1994); G. Kresse and 
D. Joubert, Phys. Rev. B 59, 1758 (1999). 
J. Iniguez and D. Vanderbilt, in preparation. 



